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The influence of a small relative density difference (~ 3.10^*) on the displacement of two miscible 
liquids is studied experimentally in transparent 2D networks of micro channels with a mean width a 
held vertically. Maps of the local relative concentration are obtained by an optical light absorption 
technique. Both stable displacements in which the denser fluid enters at the bottom of the cell and 
displaces the lighter one and unstable displacements in which the lighter fluid is injected at the 
bottom and displaces the denser one are realized. Except at the lowest mean flow velocity U, the 
average C{x,t) of the relative concentration satisfies a convection-dispersion equation. The relative 
magnitude of |?7| and of the velocity Ug of buoyancy driven fluid motions is characterized by the 
gravity number Ng = Ug/\U\. At low gravity numbers lA^gl < 0.01 (or equivalently high Peclet 
numbers Pe = Ua/Dm > 500), the dispersivities Id in the stable and unstable conflgurations are 
similar with Id oc Pe"'^. At low velocities (|A'^g| > 0.01), Id increases like 1/Pe in the unstable 
configuration {Ng < 0) while it becomes constant and close to the length of individual channels 
in the stable case (Ng > 0). Iso concentration lines c(x,y,t) — 0.5 are globally flat in the stable 
configuration while, in the unstable case, they display spikes and troughs with an rms amplitude o"/ 
parallel to the flow. For Ng > —0.2 a/ increases initially with the distance and reaches a constant 
limit while it keeps increasing for Ng < —0.2. A model taking into account buoyancy forces driving 
the instability and the transverse exchange of tracer between rising flngers and the surrounding fluid 
is suggested and its applicability to previous results obtained in 3D media is discussed. 
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I. INTRODUCTION 



Miscible displacements in porous media are encoun- 
tered in many environmental, water supply and indus- 
trial problems [H, 0, [I]- Specific types of miscible dis- 
placements, like tracer dispersion, are also usable as di- 
agnostic tools to investigate porous media heterogeneities 
at the laboratory [ij or field scales. The characteristics 
of these processes, such as the width and the geometry 
of the displacement front, are often influenced by con- 
trasts between the properties of the displacing and dis- 
placed fluids such as their density [1, H, 0, Q . In unstable 
density contrast configurations (fluid density increasing 
with height), gravity driven instabilities may appear and 
broaden the displacement front: as an unwanted result, 
this may lead to early breakthroughs of the displacing 
fluid. An example of such effects is the infiltration of a 
dense plume of pollutant into a saturated medium. 

The objective of the present paper is to study exper- 
imentally at both the local and global scales miscible 
displacements of two fluids of slightly different densities 
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{Ap/p ~ 3 X 10^^): of particular interest is the influence 
of buoyancy driven flow perturbations on the structure 
and development of the mixing zone. 

In porous media, the characteristic velocity (counted 
positively for upwards flow) of buoyancy driven flow com- 
ponents is: 



Apg 



(1) 



in which Ap is the density of the lower fluid minus that 
of the upper one, p their viscosity and k the permeability 
of the medium. The velocity Ug is the estimation, from 
Darcy's law, of the flow per unit area induced by the 
difference of the hydrostatic pressure gradients in the two 
fluids. The relative magnitude of Ug and of the mean flow 
velocity J7 is a key element in the present problem; it will 
be characterized by the gravity number [8[ : 
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With the above definition of Ap, one has Ng > in the 
stable configuration (denser fluid below the lighter one) 
and Ng < in the unstable one. 

Recent experiments 0, |6|, 0| have shown that, even 
when the parameter Ng is quite small, the geometry 
of the mixing fronts may still be strongly influenced by 
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buoyancy. Variable density flow and transport in porous 
media has therefore received an increasing attention, in 
particular through theoretical and numerical investiga- 
tions 

In the present work, optical measurements of misci- 
ble displacements of fluids of slightly different densities 
are performed in a transparent two-dimensional vertical 
network of channels with random widths [T^ . The ex- 
periments combine visualizations at the pore scale (one 
of the fluids is dyed) and measurements of the global con- 
centration profiles parallel to the mean flow at different 
mean velocities U. Comparing displacement processes 
in stable and unstable density contrast configurations al- 
lows one to detect and characterize the development and 
influence of buoyancy driven flows. 

Here, the density contrasts are low enough so that the 
\Ng\ remains small; under such conditions the hydrody- 
namic dispersion damps the development of instabilities. 
Then the mixing process can be considered as disper- 
sive and is well described by the macroscopic convection- 
dispersion equation classically used for passive tracers 
with: 

— =V.{U.C-D.VC), (3) 

where C is the tracer concentration, U the flow velocity 
and D the dispersion tensor (all values are averaged over 
the gap of the cell) ; D is assumed to reduce to the diago- 
nal components and -Dj_ corresponding to directions 
respectively parallel and perpendicular to the mean flow. 
The values of Z^n and D± are determined by two main 
physical mechanisms: advection by the velocity field in- 
side the medium and molecular diffusion (characterized 
by a molecular diffusion coefficient Dm)- The relative 
magnitude of these two effects is characterized by the 
Peclet number Pe = Ua/D„i (a is here the average chan- 
nel width). 

Both immiscible displacements [l2'| and tracer disper- 
sion [J, il3i] have already been measured previously in 
such models. This latter work [l^ uses the same ex- 
perimental technique and porous model as the present 
one but deals with experimental conditions in which the 
development of buoyancy driven instabilities is negligi- 
ble. In this case, the dye can be considered as an "ideal" 
tracer that does not modify the fluid properties. In con- 
trast, the present work deals with the influence of buoy- 
ancy effects on dispersion: the components of D depend 
on A/9 and are larger in the unstable configuratio. Similar 
studies might be performed on ZD porous samples u sing 
NMR-Imaging, CAT-Scan, acoustical techniques [ij, [l5| 
or Positron Emission Projection Imaging [l6[ but at a 
higher cost and/or with strong constraints on the fiuid 
pairs to be used. 

In the present displacement experiments, concentra- 
tion maps obtained for a vertical fiow are compared for 
different flow velocities and fluid rheologies. At the global 
scale, an effective dispersion coefficient is determined and 
its dependence on the flow velocity is studied in both sta- 



ble and unstable density contrast configurations. At the 
local scale, the variation of geometrical front features of 
different sizes is analyzed as a function of the flow ve- 
locity and of time. The combination of these local and 
global data provides both a sensitive detection of the in- 
stabilities and information on the characteristics of the 
displacement process at different length scales. The 2D 
models of porous media used here are characterized by a 
random spatial distribution of the local permeability: the 
development of the front geometry under the combined 
action of these permeability variations and of destabiliz- 
ing density contrasts is of particular interest. 

II. DESCRIPTION OF EXPERIMENT 
A. Experimental setup and procedure 

The experimental system and the technique for an- 
alyzing the data have already been described in refer- 
ence [l3|. The model medium is a vertical transparent 
two dimensional square network of channels of random 
aperture : it has a mesh size equal to c? = 1 mm and 
contains 140 x 140 channels with a mean length 0.67 mm 
and a depth 0.5 mm. The width of the channels takes 7 
values between 0.1 and 0.6 mm with a log-normal distri- 
bution and a mean value a = 0.33 mm. The permeability 
of the network is fc = 3.10"^ m^ {i.e. 3000 Darcy). 

The model is vertical with its open sides horizontal 
(see Figure 2 in reference [3). Its upper side is con- 
nected to a syringe pump sucking the fluids upward out 
of the model from a reservoir inside which the lower side 
is dipped. Initially, the model is saturated by pumping 
the first fluid of density pi out of the lower reservoir into 
the model. Then, the pump is switched off and the lower 
side of the model is removed from the liquid bath by 
lowering the reservoir (the connection tubes are shut to 
avoid unwanted fluid exchange between the model and 
the outside during this process). The reservoir is then 
emptied completely, filled up by the second fluid of den- 
sity P2 = Pi + Ap and raised again until the lower side of 
the model is below the liquid surface. The displacement 
process is initiated by opening the connection tubes and 
switching on the pump. This procedure provides a per- 
fectly straight initial front between the two fluids at the 
beginning of the displacement. 

In this work, the absolute value |Ap| of the density 
difference between the two fluids is constant with |Ap| — 
0.3 kg/m'^ . For each flow rate and pair of fluids used, both 
stable (Ap > 0) and unstable (Ap < 0) configurations 
are studied by swapping fiuids 1 and 2. 

B. Fluid characteristics 

Newtonian water-glycerol mixtures or shear thinning 
water-polymer solutions are used in the experiments. 
The mean flow velocity ranges from 0.005 to 2.5 mm.s~^. 
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The injected and displaced fluids are identical but for 
Water Blue dye added to one of the solutions, allowing 
one to both measure the local concentration optically and 
to introduce a controllable density difference between the 
fluids (note that, since the density contrast is purely due 
to the dye, the optical determination of the local dye 
concentration also measures the local density which is 
proportional to this concentration). The molecular dif- 
fusion coefficient _D,„ = 6.5 x 10~'*mm^.s~^ of the dye is 
practically the same in the polymer solutions as in pure 
water. 

The Newtonian water-glycerol solution used is ob- 
tained by mixing 60% in weight of glycerol in pure water. 
Its viscosity is equal to /i = 10^^ Pa.s at 20 °C so that 
the molecular diffusion coefficient (proportional to /i^^) 
is D*„ = 0.7 X lO-^mm^.s-^ 

The shear thinning solutions have concentrations Cp = 
500 and 1000 ppm of high molecular weight Scleroglucan 
(Sanofi Bioindustries) in water. The variation of their 
effective viscosity fi with the shear rate 7 (see Figure 1 
of ref. ^]) are well adjusted by the Carrcau function: 

^^=T-. (Mo - /^oo) +Moo. (4) 

(l + (7/7o)^) ^ 

The values of these rheological parameters for the solu- 
tions used in the present work are listed in Table [H The 
limiting value /ioo cannot be measured directly and is 
taken equal to the viscosity of the solvent, i.e. water, 
namely: fioo = 10"'^ Pa.s. In Eq. Q, T'o corresponds to 
the transition between a "Newtonian plateau" domain 
at low shear rates (7 < T'o) in which /i = /ig and a do- 
main in which fi decreases with 7 following the power law 
fi oc 7("-i) (7 > 70) ■ 

Cp n 70 Ho U* 

ppm s"'^ mPa.s mm.s~^ 

500 0.38 ± 0.04 0.077 ± 0.018 410 ± 33 0.045 
1000 0.26 ± 0.02 0.026 ± 0.004 4500 ± 340 0.013 

TABLE I; Rheological parameters of scleroglucan solutions 
used in the flow experiments. Cp refers to the polymer con- 
centration in ppm. The transition mean velocity U* is esti- 
mated by assuming a velocity gradient equal to 70 at the walls 
of a channel of width a with a parabolic velocity profile. 



C. Characteristic parameters of buoyancy driven 
flows 

For the water-glycerol mixture, all parameters are fixed 
and the modulus \Ug\ of the buoyant flow velocity (see 
Eq.(IT])) is constant (only its sign changes when the fluids 
are swapped): the gravity number Ng varies then as 
with the velocity and the influence of gravity is largest 
at the lowest velocities. 

For polymer solutions, the effective viscosity varies 
with the shear rate following Eq. (jH). and the value of 
Ng can only be computed at low ffow velocities U < U* 
for which the shear rate 7 is lower than in the whole 
model and fi = iiq — est. so that Ng varies again as 
C/~^. At velocities U > U* such that the shear rate 
is larger than 7*0 in some parts of the model, the viscos- 
ity /i varies spatially and Ng cannot be computed simply. 
However, at high shear rates one has: fi cc dotj'^~^ . Since 
Ng (X l/{fi\U\), its effective value will be proportional to 
7~": although its actual value cannot be computed, Ng 
also decreases when U increases in this shear thinning 
regime. Due to this globally monotonous variation, the 
first instabilities will still appear at the lowest velocities 
and, often, in the Newtonian domain where = /ig so 
that Ug — —Apgk/ fj,o. 

The velocity Ug has been estimated for the fluids used 
in the present work, leading to \Ug\ = 10~^, 2.10"^ and 
2.10"^ mm. s^^ for respectively the glycerol-water mix- 
ture and the 500 ppm and 1000 ppm polymer solutions 
(assuming /i = /io). At the lowest experimental flow ve- 
locity {U = 5.10^'^ mm.s^^), the duration of the complete 
saturation of the network is i ~ 3.10* s. During this time 
lapse, the distance Ug.t characterizing the growth of the 
gravitational instabilities is 28 mm for the water glycerol 
mixture: this is far above the length d of the individ- 
ual channels of the network and a noticeable influence of 
gravity on the structure of the displacement fronts is thus 
expected. For the polymer solutions, this same distance 
is less than 0.6mm, i.e. smaller than d; the dye should 
behave therefore in this case like an "ideal" tracer. 

For water-glycerol solutions, buoyancy effects should 
be sizable when the distance Ug t becomes larger than d: 
the transition should then take place at an imposed flow 
velocity Uc — Ug.L/d ~ 0.14 mm.s~^ {L is the network 
length) . This velocity correspond to the following gravity 
and Pcclct numbers : N^ ~ 0.01 an d Pe" ~ 500. These 
predictions will be shown In Section IIII CI to correspond 
well to the experimental results. 



The model is illuminated from the back by a fluores- 
cent light panel and images are acquired by a 12 bits high 
stability digital camera with a 1030 x 1300 pixels reso- 
lution (pixel size = 0.16 mm). Typically, 100 images are 
recorded for each experimental run at time intervals be- 
tween 2.5 and 700 s. The images are then translated into 
maps of the relative concentration C{x,y,t) of the two 
fluids in the model using a calibration procedure already 
described in refs. [l^, [l3| and commonly used 



III. EXPERIMENTAL RESULTS 

A. Qualitative observations of miscible 
displacements 

Figure [T] displays concentration distributions ob- 
served during displacement experiments using the water- 
glycerol mixture. In the stable configuration of Fig- 
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FIG. 1: Relative concentration maps for experiments using 
water-glycerol solutions of different densities at three differ- 



ent flow rates: (a, d) U = 0.005 mm.s 



0.2; (b, e) 



U = 0.025 mm.s^^ - l^sl = 0.04 and (c,f): U = f.25mm.s- 
- \Ng\ = S.fO"*. Experiments correspond to stable (d, e, f) 
and unstable (a, b, c) density contrast configurations at a time 
when the injected fluid occupies roughly half of the model. In 
these figures and in the following ones, darker shades corre- 
spond to the pure injected or displaced fluid and the lighter 
shade to a mixture of the two. Fluid flows are upwards with 
g pointing downwards. The field of view of the pictures is 
I53mm x f40mm. 



ures [T] d,e and f, the mean global front shape remains 
flat at all flow velocities. The overall width of the mix- 
ing zone increases however with the flow rate due to the 
development of fine structures parallel to the mean flow, 
particularly at the highest velocity (Fig. [Tf ) . 

In the unstable configuration, and at the lowest veloc- 
ity (Fig. [T^), large instability fingers with a width of the 
order of 10 to 15 mesh sizes appear and grow up to a 
length equal to that of the experimental model. For a 
velocity four times higher, fingers still appear but they 
are significantly shorter (Fig. [T|d). As the velocity in- 
creases, the size of the fingers parallel to the mean flow 
decreases while finer features develop. For a velocity still 
50 times higher (Fig. [TJ;) , the front geometry is more sim- 
ilar to that observed in the stable case although its width 
parallel to the flow is still broader. In section IIIIDl the 
structure of the mixing zone in these experiments will be 
analyzed again in the unstable case from the geometry of 
the iso-concentration fronts {C{x,y,t) = 0.5). 

Similar pictures from experiments realized with a Cp — 
1000 ppm water-scleroglucan solution are displayed in 
Figure [H In contrast with the case of the water-glycerol 
solutions, no buoyancy induced instability appears in the 



unstable configurations (cases a,b) , even at low velocities 
(a) and the relative concentration maps in both configu- 
rations are very similar. At high velocities, fine structures 
similar to those observed for the Newtonian solution ap- 
pear on the front (cases a,c). Finally, comparisons of 
front geometries observed for Cp = 500 and 1000 ppm 
and displayed in Figure 3 of ref. show that interme- 
diate scale structures of the front (with a width of 10-20 
mesh sizes) appear at the same locations for both con- 
centrations. 
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FIG. 2: Relative concentration maps for experiments using 
1000 ppm water-scleroglucan solution in gravitationally un- 
stable (a,b) resp. stable (c,d) configurations for several flow 
velocities and gravity number values ; U = 0.005 mm.s~^ - 



\Na 



4.10 (a,c) and U = 0.5 mm.s 



4.10" 



(b,d). Grey scale codes and flow direction are the same as in 
Figure [1] and the field of view is also identical. 

Globally, the above results are in qualitative agreement 
with macroscopic dispersion rneasurements realized on 
three-dimensional bead packs [5j] using conductivity trac- 
ers detected at the outlet of the samples. In these exper- 
iments, buoyancy driven instabilities are also observed at 
low velocities for Newtonian water-glycerol solutions in 
a gravitationally unstable configuration but do not ap- 
pear for water-scleroglucan solutions. Compared to this 
latter work, the optical measurements provide additional 
information on the dependence of the front structures of 
different sizes on the control parameters. We examine 
now the variation of the global dispersion characteristics 
as a function of the experimental parameters and of the 
configuration of the fluids. 



B. Quantitative concentration variation analysis 

The procedure for determining a global dispersion coef- 
ficient from the concentration maps is described in detail 
in reference p^ . The mean relative concentration C{x,t) 
of heavy (dyed) fluid at a distance x from the inlet side 
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FIG. 3: Normalized mean concentration variation C{x,t) 
as a function of the time t for unstable displacement ex- 
periments using a water glycerol mixture (C: average of lo- 
cal pixel concentration c{x,y,t) over a width Ay ~ 144 mm 
in the central part of the model), (a) U — 0.025 mm. s~^, 
Ng = -0.04, X = 104 mm. - (b) [/ = 0.005 mm.s~\ 
Ng = —0.2, X = 88 mm. Continuous line: fit by a solution 
of Equation © with (a) t = 4350 s, -D||/(7^ = 110 s and (b) 
t = 18266 s, Dii/U^ = 1369 s 

. The determination of the mean velocity U is discussed 
below in the present section IIIIBI 
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FIG. 4: Variation of the fitting parameter Z)||/(7^ as a func- 
tion of the distance x from the inlet for two unstable dis- 
placement experiments using the water glycerol solutions with 
U = 0.025 mm.s-i {Ng = -0.04) and U = 0.005 mm.s"^ 
{Ng — —0.2). Insets: variation of the mean transit time I as a 
function of x. Solid line: linear regression of the data. Values 
oft and D\\/U'^ are obtained by fitting the mean concentra- 
tion variation C{x,t) using Eq. ((3} at each distance x. 



is first determined by averaging the value c{x,y,t) for 
individual pixels over a window of width Ay = 144 mm 
representing almost the full width of the model across 
the flow (only pixels located inside the pore volume are 
included in the average). More local information is ob- 
tained from the geometry of iso-concentration curves and 
will be discussed below. 

Figures[3^-b display the variations with time of C{x, t) 
for two different values of the gravity number Ng. These 
variations have been fitted (continuous line) by the fol- 
lowing one dimensional-solution of the convection diffu- 
sion equation ^ (assuming a concentration C constant 
with y, an initial step-like variation of C at the inlet and 
taking D\\ equal to the diagonal component of the tensor 
D in the x direction) : 



C{x,t) = 



1 



1 - 



.erf 



t - t 



(5) 



Since flow is always upwards, adding the factor 
Ng/\Ng\ — ±1 makes the equation usable both for 
light fluid displacing heavy fluid (unstable configuration 



- Ng < 0) and heavy fluid displacing light fluid (stable 
configuration - Ng > 0). 

For Ng — —0.04, one observes a good fit of the exper- 
imental data with Eq. ([5]). For iV^ = —0.2 for which the 
distortions of the front due to the instabilities are very 
large, the concentration does not decrease smoothly but 
the experimental curve displays bumps; these features 
are the signature of rising fingers reaching the measure- 
ment height. Yet, an acceptable fit of the experimental 
data with Eq. ([5]) can still be achieved. The fits provide 
the values of the mean transit time t of the front at the 
distance x and of the ratio D\\/U^ — At^/(2Z) (Ai^ is 
the centered second moment of the transit times along 
the distance x). 

The insets of figures [5^-b display the variation of t 
with the distance x which is linear in both cases; this 
shows that the mixing zone characterized by the mean 
concentration profile C'{x, t) moves at a constant velocity 
Umz which is equal to the inverse of the slope of the 
variation. A linear regression of the data provides the 
values Umz — 0.027 mm.s~^ for case (a) and Umz = 
0.0046 mm.s^^ for (6). The mean fiuid velocity U in the 
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model may be taken equal to the ratio of the injected 
flow rate and of the pore volume per unit length along x 
which have been determined independently. The values 
of U computed in this way at the same two flow rates as 
above are respectively U — 0.025 and 0.005mm.s~^: they 
are very close to the corresponding values of Umz- 

The two velocities U and Umz are compared more pre- 
cisely in Figure [S] which displays the values of Umz/U 
for both stable (o) and unstable (•) flow configurations,. 
This ratio is always close to 1: this shows that buoy- 
ancy effects do not influence the mean displacement of 
the concentration profile, even at the lowest flow rate 
{Ng = —0.2) for which they are very large (Figl2^). One 
assumes therefore in the following that U = Umz- 
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FIG. 5: Variation as a function of |A^g| of the characteristic 
velocities of the fluid displacement normalized by the mean 
velocity U in stable (open symbols) and unstable (dark sym- 
bols) density contrast configurations. (o),(»): normalized ve- 
locity of the mixing region Umz/U determined by a linear re- 
gression on the variation of I as a function of x (see insets of 
Fig.HK-b). (□, B) UcO.s, ~ mean velocity of iso-concentration 
lines c — 0.5. 

The main graphics of Figures [4l^a-b) display the vari- 
ation of the ratio D\\/U'^ as function of the distance x 
from the inlet. In case (a) {Ng = -0.04), £'||/C/^ in- 
creases at first slightly with the distance x and levels off 
for a; > 50 mm; then, it fluctuates around a constant 
value. Previous studies [l^ have shown that the fluctua- 
tions are periodic and determined by the structure of the 
network (the period is equal to the mesh size). In case (b) 
{Ng — —0.2), fluctuates around a constant value 

for X < 100 mm. At larger distances, a slight steady in- 
crease of D\\/U'^ with X is observed. This is likely due to 
the rising fingers, directly observable in Fig[T^ and which 
can be identified on the curve of FiglSb. Tracers are ad- 
vected faster and farther inside these fingers than outside, 
resulting in an increase of the dispersivity. Even in this 
case, however, the fluctuations of the measurements and 
the relatively small variations of the values of D\\/U'^ do 
not allow one to conclude that a diffusive regime is not 
reached. 



These results show both that the mixing front moves 
at a constant velocity U and that (except perhaps for 
Ng — —0.2) it reaches a dispersive spreading regime char- 
acterized by a dispersion coefficient < D[| > taken equal 
to the average of over the full experimental range of x 
values. In the following, like in ref. [l3l |. the dispersivity 
Id = < D\\ > /U is generally used instead of < Dy > to 
characterize the dispersion process. Finally, the standard 
deviation of the individual values of Z3|| will be used to 
estimate the error bars on l^- 



C. Global dispersion measurement results 
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FIG. 6: Variation of the dispersivity /^/a = D/aU, as a 
function of the Peclet number Pe for different fluids, (o, •): 
water- glycerol mixture; (□, ■): 500 ppm and (A, A): 1000 
ppm polymer solutions. Light symbols: Ng < 0; dark sym- 
bols: Ng > 0. Upper axis: values of |A'^g| for the New- 
tonian water-glycerol solutions. Dashed line: variation law 
Id/a, — 1.5 Pe" "*^ satisfled from ref. [ij by data corresponding 
to the polymer solution. Solid {Id/a = 0.6 Pe"'^) and dash- 
dotted {Id/a = 0.5 Pe°-^) lines:, regression for water-glycerol 
data respectively for the unstable and stable configurations 
in the buoyancy free flow domain (|A^g| < Ng = 0.01). Dot- 
ted line: variation predicted by Eq. (|10|l based on the model 
described in section ITVl below (the dispersivity is assumed to 
be the sum of a passive tracer component and of a buoyancy 
term which decreasing with Pe). The fitting parameters are 
/3t = and e = 0.3. 

In the present section, we discuss the variations of the 
normalized dispersivity Id/o, as a function of the flow ve- 
locity in the stable and unstable configurations: these 
variations are displayed in Figure [S] for the different solu- 
tions investigated (the horizontal scales are the Peclet 
number Pe = Ua/Dm (bottom axis) and the gravity 
number \Ng\ (top axis). 

For water-glycerol solutions, there is a clear separation 
at high Ng (or equivalently low Pe) values between data 
points corresponding to unstable {Ng < 0) and stable 
{Ng > 0) density contrast configurations: this separation 
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occurs occurs close to the transition value = 0.01 
(respectively Pe'^ = 500) discussed in Sec. Ill CI Such 
a separation has also been previously reported in 3D 
bead packings 5]: in Sec. IIVI these 31? observations will 
be compared quantitatively to the present 2D measure- 
ments. This separation reflects the development of finger 
like structures at low velocities in the unstable configura- 
tion {Ng < 0) and the flattening of the front in the stable 
one {Ng > 0). These features are due to buoyancy forces 
and are qualitatively visible in Figure O At the lowest 
Peclet number, the values of Id for Ng > and Ng < 
differ by a factor of nearly 10. 

For Ng > 0, the dispersivity Id reaches for Ng ~ 2 x 
10^'^ [Pe ~ 100) a low minimum value Id — ^ mm close 
to the characteristic local length, i.e. the mesh size of the 
lattice. In most usual homogeneous porous media, and 
in the case of a perfectly passive tracer, this latter value 
represents a lower limit of Id', finding such a low value 
here confirms the stabilizing influence of the buoyancy 
forces ^8]. 

For displacements of water-polymer solutions in the 
same range of Pe values, the variations of Id with Pe 
are the same for the stable and unstable flow conflgura- 
tions: this result was qualitatively visible in Figure [2] for 
the 1000 ppm solution and is confirmed quantitatively in 
Figure [6] for both polymer concentrations. This implies 
that no fingering instabilities develop in the unstable case 
and that the front is not fiattened by gravity in the stable 
one. 

Quantitatively, from Sec. Ill CI an upper limit of the 
value of Ng can be estimated by taking ^ = /xq in Eq. ^ 
at the lowest experimental flow velocity. This gives: 



5 X 10 and |7V„| = 4.5 x 10 respectively for 



the 500 and 1000 ppm solutions, i.e. below the thresh- 
old of the instabilities observed for the Newtonian fluids. 
Since, Ng decreases at higher flow velocities, this con- 
firms that no effect of buoyancy will be observable in 
the experimental range of velocities (in agreement with 
the experimental observations and with the discussion of 
SecHTCl). 

At high velocities {U > 0.1 mm. s^^ (or Pe > 500 for 
water-glycerol and Pe > 50 for the polymer solutions), 
Id takes similar values for both Ng > and < and 
increases with Pe as Id oc Pe" with a — 0.5 at high 
Pe values in both the unstable and stable configurations 
(solid and dash dotted lines in Figure E). A similar be- 
haviour has aleady been reported ref. [13 1 for water poly- 
mer solutions but with a lower exponent a — 0.35: the 
variation differs in both cases from the slow increase ob- 
served in three-dimensional media such as an homoge- 
neous monodisperse grain packing Q . The Pe" variation 
reflects the combined effects of geometrical dispersion due 
to the disorder of the velocity field and of Taylor disper- 
sion due to the velocity profiles in individual channels. 
The latter becomes important at high Peclet numbers 
due to the reduced mixing at the junctions by transverse 
molecular diffusion. As a result, the correlation length 
of the velocity of the tracer particles along their trajec- 



tories (and, therefore, the dispersivity Id) become larger. 
The difference between the values of a for the Newtonian 
and shear thinning solutions is likely due to the different 
relative weight of the two mechanisms in the two cases: 
as discussed in refs. [H, Taylor dispersion is indeed 
reduced and geometrical dispersion is enhanced for shear 
thinning fluids compared to the Newtonian case. 

The above analysis of the global dispersion discussed 
above uses averages of the local concentration over nearly 
the full width of the model and encompassing therefore 
all the geometrical features of the front instabilities. The 
corresponding global value of Id combines then different 
types of effects. 

A first one is the local spreading of the displacement 
front. It is likely to result from the combined effects 
the local disorder of the flow fleld (geometrical dispersion 
mechanism) and the flow profile between the rough walls 
(Taylor mechanisms). The effect of these mechanisms in 
the present system is discussed in ref. [l^ 

A second, more global, effect is the global spreading 
of the mixing zone due to the fluid velocity contrasts 
between different flow paths (associated, for instance, to 
the development of the instabilities) . This is studied now 
specifically from the variations of the front geometry. 



D. Spatial structure of the displacement fronts for 
unstable flows 



Previous experiments in three dimensional porous me- 
dia ,15i demonstrated a clear amplification of front struc- 
tures resulting from permeability heterogeneities for un- 
stable density contrast configurations. Such effects can 
be studied precisely here down to small lengths scales 
thanks to the two-dimensional geometry of the model 
network and to the high precision and spatial resolution 
of the optical concentration measurements. 

In the following, the front geometry is characterized 
from the lines Xf{y, t) along which the local relative con- 
centration c(x, y, t) is equal to 0.5 at a given time t. Ex- 
amples of such lines determined by a thresholding proce- 
dure at four gravity numbers —0.001 > Ng > —0.2 are 
displayed in Fig. \7\ 



At the lowest flow velocity investigated {Ng 



-0.2, 



Pe ~ 25), the buoyancy driven flow components have 
a major influence on the front geometry and several in- 
stability flngers soar up while the front displacement is 
much slower in other regions (Fig.[7li). 

At higher mean flow velocities (Fig. EId-c) the front re- 
tains a rough geometry, still reflecting buoyancy driven 
flow components but its mean advancing motion is clearly 
visible. In contrast, for Ng = —0.2 (Fig. [TJi), the devel- 
opment of the front appears as the combination of the 
independent growth of the individual flngers. 

Another important feature is the fact that the distance 
y across the front at which corresponding geometrical 
features (peaks and troughs) appear remains the same. 



For No 



-0.2, the extension of these features parallel 



8 



N =0.00125 



"1 — ' — r 



N =0.02 



50 lOO 150 



50 100 150 



150 



N =0.04 



J I L 



1 


1 1 1 

N,=0.2 


















1 







50 lOO 150 50 lOO 150 

FIG. 7: Iso concentration fronts Xf{y,t) measured at four 
different gravity numbers Ng = -0.00125 (a), -0.02 (b), 
—0.04 (c) and —0.2 (d) for water-glycerol mixtures in an un- 
stable flow configuration (the injected fluid occupies half of 
the model area). Flow is upwards with g oriented downwards. 
The inset in figure b, displays the typical size of the spikes A. 



to the flow increases and they cluster together into larger 
structures as can be seen in Figs. [71d-c. Moreover, for 
-0.02 > Ng > -0.2, and even though the width of the 
front parallel to the mean flow increases significantly, the 
typical transverse size A (along y) of the individual spikes 
of the front is fairly constant wih A ~ 4 — Smm. This 
value has been taken equal to the mean interval (along 
y) between successive local maxima of the local distance 
(parallel to x) of the front from the inlet side (see insert 
on figure Hb). 

These results contrast with the assumptions of a vary- 
ing wavelength P3| often applied to porous media follow- 
ing observations in Hele-Shaw cells ^20]. In the Hele-Shaw 
case, there is however no characteristic length scale of the 
geometry of the flow field inside the cell, beyond its thick- 
ness; in the present case, on the contrary, the location of 
the features of the front appears therefore be determined 
by the heterogeneities of the flow field. 

At still higher flow velocities (for instance in Fig. [7^) 
and for Ng > -N;; = -0.01 (See Sec. IlC] for the ex- 
pression of Ng), distortions of the front due to buoyancy 
effects decrease in size and become hardly visible on the 
iso-concentration lines. 

Quantitatively, we characterize these iso-concentration 



fronts by their mean position along the flow Xf{t) =< 
Xf{y, t) >y and by the rms fluctuations cry of the distance 

As shown in the inset of Figure [51 the mean distance 
Xf{t) increases linearly with time even for the lowest 
flow velocity Ng — —0.2. The propagation of the iso- 
concentration fronts may therefore be characterized by a 
velocity Uco.5 determined from a linear regression on the 
variation of Xf{t) with t (□ and ■ symbols in Fig. O. 
The velocity C/cO.5 is very close to the flow velocity U 
(see Sec HirCI) for Ng > -0.08; however, it is 40% higher 
at the lowest flow velocity {Ng = —0.2). This confirms 
other indications of a transition towards a different type 
of front growth dynamics (in this latter case, the value of 
Uco.5 is likely determined by the development of the few 
large fingers observed in Figure [8]). 
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FIG. 8: Variation of the fronts width af as function of the 
mean front position xj. Empty symbols stand for stable 
configurations, filled symbols for unstable ones. Squares, 
triangles, diamonds and circles correspond respectively to 
0.04 and 0.2. Inset: variation of xJ as a 



function of time in the unstable configuration for Ng 



-0.2 



(•), Ng = -0.08 (*) and iV^ 
regressions on data points. 



-0.04 (♦). Solid lines: linear 



The variation of cry with Xf{t) depends also strongly 
on the value of Ng (Fig. [5]). At the highest mean flow 
velocity {i.e. \Ng\ = 0.00125), ct/ becomes constant and 
equal to 6 mm as soon as the distance from the injection 
line is larger than 10 mm. This value is the same for 
both stable and unstable flow configurations: this is in 
agreement with the qualitative observations of Sec. IIII Al 
in which no influence of buoyancy is visible when com- 
paring Figs.[TJ; and [If. 

In contrast, at low flow rates (or large |iVg| values), 
the large differences between the corresponding concen- 
tration maps in Fig. [T] are reflected in the variations of ct/ 
with x/. At the lowest flow-rate {\Ng \ = 0.2), a/ keeps in- 
creasing roughly linearly with xJ in the unstable case; in 
the stable configuration, in contrast, cr/ reaches quickly a 
limit of the same order of magnitude as for the highest ve- 
locity. At the intermediate velocities {Ng = 0.02 — 0.04), 
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the limiting value of a/ is higher in the unstable config- 
uration than both in the stable one and at the highest 
velocity. The distance xj required to reach this limit is 
also increased. 

This result suggests that, in the unstable configuration 
and for intermediate Ng values (—0.01 > Ng > —0.1), 
two distinct regimes are successively observed. At short 
distances, the displacement is controlled by the instabil- 
ity while the variation of cr^ is similar to that measured 
for Ng — —0.2. At larger distances, the variation of ct/ 
levels off and it reaches a nearly constant value like for 
stable displacements. The transition distance increases 
with Ng-. for Ng = —0.2, it is of the order of the sam- 
ple length, explaining why the second regime is not ob- 
served. Together with results reported previously, we use 
biw this observation to provide an analytical estimation 
of the dispersivity for Ng < —Ng. 

In order to compare the results of the present sec- 
tion IIII Dl with those of section IIII C[ it must be empha- 
sized that (Tf refers to the extension of the iso concen- 
tration front c = 0.5 in the fiow direction: it docs not 
include therefore the infiuence of the width of the con- 
centration profile at a given transverse distance y. At 
high velocities, for instance, the iso concentration front 
is nearly fiat (Fig. [7^) and displays only a few spikes. In 
this case, cr/ reaches a non-zero limit at long distances 
(due to these small features) and it does not increase as 
like the global width of the mean concentration pro- 
file discussed in Sec lIII Cl In the unstable cases and at 
low velocities, the linear increase of cr/ at long distances 
reflects directly the buoyant rise of fingers: its dynamics 
differs from that of the global spreading of the mixing 
zone. The latter results from the combination of several 
mechanisms (including, but not exclusively, the growth 
of the fingers) leading to an increase of the global width 
of the concenration profile as 



IV. ESTIMATION OF DISPERSIVITY 
VARIATIONS FOR UNSTABLE 
DISPLACEMENTS 

The development of fingers driven by buoyancy forces 
in the unstable configuration is opposed by lateral mix- 
ing induced by transverse dispersion: it reduces the local 
density contrast Sp between the fingers and the surround- 
ing fiuid and, finally, the buoyancy forces. As pointed out 
in Sec lIII Dl the structure of the displacement front has 
a characteristic size A constant and close to 4 mm for 
gravity numbers in the range —0.01 > Ng > —0.2. The 
characteristic exchange time r for lateral mixing will then 
be of the order of the transverse diffusion time across the 
half- width A/2 with: 



A^ 



8Di 



(6) 



The rising motion of the fingers driven by buoyancy 
forces should have a velocity u finger proportional to the 
local density contrast with u finger {^pI ^pUg)- As- 
suming that 5p decreases exponentially with a time con- 
stant of the order of the transverse mixing time t leads to 
Ufinger U g exp{—t / t) . The vertical displacement 1(t) 
of the rising finger before its velocity goes to zero should 
then be given by: 



1{t) ^ rUg 



(7) 



1{t) is then the typical spreading distance of the front 
due to buoyancy driven motions during the time r: it 
should be of the same order of magnitude as the limiting 
value of the width cr/ of the iso concentration lines at 
long distances. The transverse dispersion coefficient D±^ 
generally decreases with the Peclet number (or equiva- 
lently with the fiow velocity U) so that both r and 1{t) 
should increase at low velocities. This agrees with the 
observed increase of cr/ at long distances. At very low 
velocities, 1{t) will reache (or exceed) the system size: 
this is indeed observed at the lowest experimental fiow 
velocity {Ng — —0.2). In that case, cr/ keeps increasing 
with distance and does not reach a constant value within 
the model length. 

In order to estimate the dispersion coefficient com- 
ponent Dbuoyancy associatcd to these buoyancy driven 
motions, one may consider r represents as the char- 
acteristic crossover time towards diffusive front spread- 
ing: the corresponding width Z(t) should then verify 
Z(t)^ ~ 2DiiuoyancyT. Combining with Eqs. ([6]) and (O, 
this leads to the following value for the dispersivity com- 
ponent buoyancy ' 



Id 



D 



buoyancy 



buoyancy 

U 



16UDj_ ■ 



(8) 



As a first approximation, the total dispersivity Id is con- 
sidered to be the sum of the dispersivity for a fully 
passive tracer and of Id buoyancy- this amounts to as- 
sume that these two dispersion processes are indepen- 
dent. This assumption us only an approximation since 
the spatial fiow velocity variations in the model infiuence 
both spreading processes. 

The normalized passive tracer dispersivity Z^'^/a may 
be estimated from data obtained with the polymer so- 
lutions for which, as noted above, no buoyancy effect is 
visible. The equation: 



eq 



fPe' 



(9) 



where a = 0.35 ± 0.03 and / = 1.5 ± 0.3 has been se- 
lected for that purpose: it provides indeed a good global 
fit (dashed line in Figure [6]) with the polymer data at all 
Pe values and an acceptable one with water-glycerol data 
either in the stable configuration or at high velocities. 
Combining Eqs. ([Hini) and replacing U by its expression 
function of Pe leads then to: 



in which D± is the transverse dispersion coefficient. 



Id 
a 



2\2 



/Pe" 



up 



16ePe^+0TDl 



(10) 
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in which e is a constant. In this equation the transverse 
dispersion coefficient D± is assumed to vary as: D± — 

□ 




: ffi 



1 1- 1 

0.01 0.02 0.03 0.04 0.05 

1/Pe 



report to be equal and to range between 1 and 1.2 [3[. 
As a result, the difference Id — l^'^ estimated from the 
model is expected to increase with Pe like ~ 1/Pe^. We 
have tested this prediction is tested by re-analyzing the 
dispersion coefficients measured by Freytes et al [5]. In 
this work, gravity driven instabilities were studied in a 
model porous packing of 1 mm glass beads and for a 
density difference Ap = 10~^g/l. For water and in the 
high Pe regime where the effect of buoyancy is negligible, 
their data show that : l'^'^ oc Pg0.i±o.i indicating that 
/3t — 1-1 ± 0.1. As above, the buoyancy component is 
estimated by subtracting the passive tracer dispersivity 
ZJ^' from the measured dispersivities. The difference — 
obtained in this way is found to vary as (1/Pe)^'^^°'^ 
(see inset of Fig. ^ with an exponent 1.9 close to the 
value 1-|-/3t — 2.1±0.1. predicted by the model 



V. CONCLUSIONS 



FIG. 9: Variation of Id/ a as function of 1/ Pe for the water- 
glycerol solution. Filled (resp. empty) circles correspond to 
Ng < (resp. Ng > 0). Inset : Variation of {Id — I'd^)/^ as 
function of 1/Pe for unstable experiments in a 2D network 
(•) and a 3Z) porous medium (■) [^. In order to make com- 
parisons easier, horizontal (resp. vertical) values for 3D data 
have been divided by a factor 5 (resp. 0.6). Solid (resp. dot- 
ted) lines: regressions using power laws of exponents 1 (resp. 
2.1). 

In order to put more emphasis on the buoyancy con- 
trolled regime at low flow velocities (low Pe), the vari- 
ation of the dispersivity is plotted in Figure [9] as func- 
tion of 1/Pe. For unstable flows, as soon as 1/Pe > 
0.0025 (1/400), Id steadily increases with 1/Pe. This is 
also the case of the buoyancy component estimated by 
subtracting the passive tracer dispersivity component Z^' 
from the values of Id- The difference Id — l*/^ increases 
linearly with 1 / Pe and can be fltted by the second term 
of Eq. pU]) with e = 0.3 and /?t — (continuous line in 
the inset of Figure [H]) . 

The experimental value of Pt is lower than that re- 
ported in ref. 21] {(3t = 0.2); this result corresponds to 
numerical simulations for capillary tube networks with a 
normalized standard deviation Ua/a of the channel aper- 
ture equal to that of ours {aa/a, = 0.3). However, in 
our experiments, the influence of the Id buoyancy term is 
mostly significant at the lower Peclet numbers (Pe < 50) 
while the simulations of ref. [2l[ deal with Pe > 300. The 
variation of D±^ for Pe < 50 may therefore be expected 
to be slower (and the corresponding exponent lower) due 
to the influence of the molecular diffusion coefficient Dm 
which represents a constant lower limit at very low Peclet 
numbers. 

In the present work, the variation of the coefficient of 
dispersion D results directly from the 2D nature of the 
network which has been used. In a 3D porous medium, 
a bead packingfor instance, D depend weakly on the 
Peclet number [3, H, [1] . In this case Pt and a are usually 



To conclude, the miscible vertical displacement mea- 
surements on transparent networks of channels reported 
here have provided information at both the local and 
macroscopic scales on the mixing front of two misci- 
ble fluids of slightly different densities. Qualitatively, 
the present macroscopic dispersion measurement confirm 
previous ones performed on three-dimensional porous 
media [8.] : the key feature of this work is that additional 
new information is provided by the high resolution visu- 
alization of front structures of different sizes down to the 
scale of individual channels. Using this information, the 
front distortions resulting from instabilities in unstable 
density contrast configurations could be analyzed quan- 
titatively. 

In these systems, the global spreading of the front re- 
sults from the combination of the effects of the disordered 
spatial variations of the velocity field (only mechanism 
active in the passive tracer case) and of buoyancy driven 
flows; the latter may either decrease or reduce the dis- 
persion depending on the gravitationally stable or unsta- 
ble conflguration of the fluids. For large Peclet numbers 
(Pe > 500) [i. e. small gravity numbers \Ng\ < 0.01), 
the displacement fronts are very similar in both config- 
urations and the global shape of the iso-concentration 
fronts is flat. In this case, the front spreading charac- 
teristics and the dispersivity values are similar to those 
measured for shear-thinning polymer solutions. 

In the stable configuration, the dispersivity decreases 
significantly for Pe < 500 and becomes constant and 
close to 1 mm when Pe < 100. At the same time, the 
geometry of the iso-concentration front is only weakly 
affected by the stabilizing effect of buoyancy. 

In the unstable configuration, and at moderate Ng val- 
ues (e.g. —0.01 > Ng > —0.2), the initial development of 
the instabilities is damped by transverse hydrodynamic 
dispersion after some time. As a result, front spreading 
remains dispersive but with a dispersivity increasing at 
low velocities as Id — 1/Pe. These instability fingers are 
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reflected in the geometry of the iso-concentration fronts 
which display large spikes and troughs. In this range 
of Ng values, the location (in the transverse direction) 
and the width of the spikes are observed to be constant. 
The spreading in the flow direction as a function of time 
displays two regimes: initially (short distances from the 
injection), the width of the iso-concentration fronts in- 
creases linearly and, then, it levels-off towards a constant 
value at large distances. 

In order to explain these results, an approach combin- 
ing the influences of longitudinal buoyancy forces parallel 
to the mean flow and of the exchange of solute between 
the instability flngers and the surrounding fluid has been 
developed. It allows one to account semi-quantitatively 
for the dependence of the dispersivity on the Peclet num- 
ber if Ng > —0.2. The present experimental observations 
on 2D networks as well as previous measurements on 3-D 
bead packings (for which Id — 1/Pe^) are well fitted by 
the model. 

At the lowest fiow velocity investigated {Ng — —0.2), 
large fingers develop on the interface; the concentration 
front is strongly distorted but its spreading can still be 
considered as dispersive. Yet, the number of fingers on 
the iso-concentration front decreases with time while its 
width parallel to the flow keeps increasing with distance. 



This suggests that, above this gravity number, grav- 
itational instabilities might control the transport pro- 
cess. In 3D porous media, a linear growth of the mix- 
ing zone reflecting the dominant influence of such insta- 
bilities is only observed above a higher threshold value 
Ng -1.5 0. 

Further studies will be needed to confirm our observa- 
tion by using pairs of fluids with different density con- 
trasts. Another issue of practical interest is the influence 
of the viscosity contrast on the spatial distribution of the 
tracer. 
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